#include "StdAfx.h"

#include <limits>
#include <math.h>
#include <vector>
#include "CalcThreshold.h"


//---------------------------------------------------------------------------
TCalcThreshold::TCalcThreshold(long* HistVal)
{
	Threshold = 0;
	for(int i=0;i<256;i++){
		this->HistVal[i] = HistVal[i];
	}
}
//---------------------------------------------------------------------------

TCalcThreshold::~TCalcThreshold(void)
{
}


//---------------------------------------------------------------------------
void  TCalcThreshold::CalculateThreshold(EnumThreshold ThresholdMethod)
{
    switch (ThresholdMethod)
    {   
       case MeanThreshold:
            Threshold = GetMeanThreshold(HistVal);
            break;

        case PTileThreshold:
            Threshold = GetPTileThreshold(HistVal);
            break;

		case MinimumThreshold:
			Threshold = GetMinimumThreshold(HistVal);
		    break;

        case IntermodesThreshold:
            Threshold = GetIntermodesThreshold(HistVal);
	        break;

        case IterativeBestThreshold:
            Threshold = GetIterativeBestThreshold(HistVal);
            break;

        case OSTUThreshold:
            Threshold = GetOSTUThreshold(HistVal);
            break;

        case OneDMaxEntropyThreshold:
            Threshold = Get1DMaxEntropyThreshold(HistVal);
            break;

        case MomentPreservingThreshold:
            Threshold = GetMomentPreservingThreshold(HistVal);
            break;

        case HuangFuzzyThreshold:
            Threshold = GetHuangFuzzyThreshold(HistVal);
            break;

        case KittlerMinError:
           Threshold = GetKittlerMinError(HistVal);
           break;

        case IsoDataThreshold:
            Threshold = GetIsoDataThreshold(HistVal);
            break;

        case ShanbhagThreshold:
            Threshold = GetShanbhagThreshold(HistVal);
            break;

        case YenThreshold:
           Threshold = GetYenThreshold(HistVal);
           break;

        case NotProcessImageBinarize:
            Threshold = 0xFF;
            break;

        default:
            break;
    }
}
//---------------------------------------------------------------------------
long TCalcThreshold::GetThreshold(EnumThreshold ThresholdMethod)
{
	CalculateThreshold(ThresholdMethod);
    return Threshold;
}
//---------------------------------------------------------------------------
long TCalcThreshold::GetMeanThreshold(long HistGram[])
{    long Sum = 0, Amount = 0;    for (int Y = 0; Y < 255; Y++){		Amount += HistGram[Y];		Sum += Y * HistGram[Y];	}	return Sum / Amount;}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetPTileThreshold(long HistGram[], int Tile)
{    long Y, Amount = 0, Sum = 0;    for (Y = 0; Y < 256; Y++) Amount += HistGram[Y];    for (Y = 0; Y < 256; Y++)    {        Sum = Sum + HistGram[Y];        if (Sum >= Amount * Tile / 100) return Y;    }    return -1;}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetMinimumThreshold(long HistGram[])
{    int Y, Iter = 0;    double HistGramC[256];    double HistGramCC[256];    for (Y = 0; Y < 256; Y++)    {        HistGramC[Y] = HistGram[Y];        HistGramCC[Y] = HistGram[Y];    }
    while (IsDimodal(HistGramCC) == false)                                         {        HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;                    for (Y = 1; Y < 255; Y++)            HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;          HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;                for(int i=0; i<256; i++)        {            HistGramC[Y] = HistGramCC[i];        }        Iter++;        if (Iter >= 1000) return -1;    }    bool Peakfound = false;    for (Y = 1; Y < 255; Y++)    {        if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peakfound = true;        if (Peakfound == true && HistGramCC[Y - 1] >= HistGramCC[Y] && HistGramCC[Y + 1] >= HistGramCC[Y])            return Y - 1;    }    return -1;}
//---------------------------------------------------------------------------
bool  TCalcThreshold::IsDimodal(double HistGram[]) 
{    int Count = 0;    for (int Y = 1; Y < 255; Y++)    {        if (HistGram[Y - 1] < HistGram[Y] && HistGram[Y + 1] < HistGram[Y])        {            Count++;            if (Count > 2) return false;        }    }    if (Count == 2)        return true;    else        return false;}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetIntermodesThreshold(long HistGram[])
{
    int Y, Iter = 0, Index;
    double HistGramC[256];       
    double HistGramCC[256];    
    for (Y = 0; Y < 256; Y++)
    {
        HistGramC[Y] = HistGram[Y];
        HistGramCC[Y] = HistGram[Y];
    }

    while (IsDimodal(HistGramCC) == false)                                                 
    {
        HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;              
        for (Y = 1; Y < 255; Y++)
            HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;     
        HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;     

        for(int i=0; i<256; i++)
        {            HistGramC[Y] = HistGramCC[i];        }  

        Iter++;
        if (Iter >= 10000) return -1;                                                    
    }

    int Peak[2];
    for (Y = 1, Index = 0; Y < 255; Y++)
        if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peak[Index++] = Y - 1;
    return ((Peak[0] + Peak[1]) / 2);
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetIterativeBestThreshold(long HistGram[])
{
    int X, Iter = 0;
    int MeanValueOne, MeanValueTwo, SumOne, SumTwo, SumIntegralOne, SumIntegralTwo;
    int MinValue, MaxValue;
    int Threshold, NewThreshold;
    for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
    for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
    if (MaxValue == MinValue) return MaxValue;       
    if (MinValue + 1 == MaxValue) return MinValue;   
    Threshold = MinValue;
    NewThreshold = (MaxValue + MinValue) >> 1;
    while (Threshold != NewThreshold)    
    {
        SumOne = 0; SumIntegralOne = 0;
        SumTwo = 0; SumIntegralTwo = 0;
        Threshold = NewThreshold;
        for (X = MinValue; X <= Threshold; X++)         
        {
            SumIntegralOne += HistGram[X] * X;
            SumOne += HistGram[X];
        }
        MeanValueOne = SumIntegralOne / SumOne;
        for (X = Threshold + 1; X <= MaxValue; X++)
        {
            SumIntegralTwo += HistGram[X] * X;
            SumTwo += HistGram[X];
        }
        MeanValueTwo = SumIntegralTwo / SumTwo;
        NewThreshold = (MeanValueOne + MeanValueTwo) >> 1;   
        Iter++;
        if (Iter >= 1000) return -1;
    }
    return Threshold;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetOSTUThreshold(long HistGram[])
{
    int X, Y, Amount = 0;
    int PixelBack = 0, PixelFore = 0, PixelIntegralBack = 0, PixelIntegralFore = 0, PixelIntegral = 0;
    double OmegaBack, OmegaFore, MicroBack, MicroFore, SigmaB, Sigma;            
    int MinValue, MaxValue;
    int Threshold = 0;
    for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
    for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
    if (MaxValue == MinValue) return MaxValue;         
    if (MinValue + 1 == MaxValue) return MinValue;     
    for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];    
    PixelIntegral = 0;
    for (Y = MinValue; Y <= MaxValue; Y++) PixelIntegral += HistGram[Y] * Y;
    SigmaB = -1;
    for (Y = MinValue; Y < MaxValue; Y++)
    {
        PixelBack = PixelBack + HistGram[Y];
        PixelFore = Amount - PixelBack;
        OmegaBack = (double)PixelBack / Amount;
        OmegaFore = (double)PixelFore / Amount;
        PixelIntegralBack += HistGram[Y] * Y;
        PixelIntegralFore = PixelIntegral - PixelIntegralBack;
        MicroBack = (double)PixelIntegralBack / PixelBack;
        MicroFore = (double)PixelIntegralFore / PixelFore;
        Sigma = OmegaBack * OmegaFore * (MicroBack - MicroFore) * (MicroBack - MicroFore);
        if (Sigma > SigmaB)
        {
            SigmaB = Sigma;
            Threshold = Y;
        }
    }
    return Threshold;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::Get1DMaxEntropyThreshold(long HistGram[])
{    int  X, Y,Amount=0;    double HistGramD[256] = {0.0};    double SumIntegral, EntropyBack, EntropyFore, MaxEntropy;    int MinValue = 255, MaxValue = 0;    int Threshold = 0;
    for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;    for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;    if (MaxValue == MinValue) return MaxValue;             if (MinValue + 1 == MaxValue) return MinValue;    
    for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];     
    for (Y = MinValue; Y <= MaxValue; Y++)   HistGramD[Y] = (double)HistGram[Y] / Amount+1e-17;
    
    MaxEntropy = DBL_MIN;    for (Y = MinValue + 1; Y < MaxValue; Y++)    {        SumIntegral = 0;        for (X = MinValue; X <= Y; X++) SumIntegral += HistGramD[X];        EntropyBack = 0;        if(SumIntegral == 0){            if(SumIntegral == 0) continue;        }        for (X = MinValue; X <= Y; X++) {            double tmp = HistGramD[X] / SumIntegral ;            if (tmp<=0) continue;            EntropyBack += (-HistGramD[X] / SumIntegral * log(tmp));        }        EntropyFore = 0;        for (X = Y + 1; X <= MaxValue; X++){            if(SumIntegral == 1) continue;            double tmp = HistGramD[X] / (1 - SumIntegral) ;            if (tmp<=0) continue;            EntropyFore += (-HistGramD[X] / (1 - SumIntegral) * log(tmp));        }        if (MaxEntropy < EntropyBack + EntropyFore)        {            Threshold = Y;            MaxEntropy = EntropyBack + EntropyFore;        }    }    return Threshold;}

//---------------------------------------------------------------------------
unsigned char   TCalcThreshold::GetMomentPreservingThreshold(long HistGram[])
{
    int X, Y, Index = 0, Amount=0;
    double Avec[256];
    double X2, X1, X0, Min;
    for (Y = 0; Y <= 255; Y++) Amount += HistGram[Y];     
    for (Y = 0; Y < 256; Y++) Avec[Y] = (double)A(HistGram, Y) / Amount;       // The threshold is chosen such that A(y,t)/A(y,n) is closest to x0.
    // The following finds x0.
    X2 = (double)(B(HistGram, 255) * C(HistGram, 255) - A(HistGram, 255) * D(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
    X1 = (double)(B(HistGram, 255) * D(HistGram, 255) - C(HistGram, 255) * C(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
    X0 = 0.5 - (B(HistGram, 255) / A(HistGram, 255) + X2 / 2) / sqrt(X2 * X2 - 4 * X1);
    for (Y = 0, Min = (std::numeric_limits<double>::max)(); Y < 256; Y++)
    {
        if (abs(Avec[Y] - X0) < Min)
        {
            Min = abs(Avec[Y] - X0);
            Index = Y;
        }
    }
    return (byte)Index;
}
//---------------------------------------------------------------------------
double  TCalcThreshold::A(long HistGram[], int Index)
{
    double Sum = 0;
    for (int Y = 0; Y <= Index; Y++)
        Sum += HistGram[Y];
    return Sum;
}
//---------------------------------------------------------------------------
double  TCalcThreshold::B(long HistGram[], int Index)
{
    double Sum = 0;
    for (int Y = 0; Y <= Index; Y++)
        Sum += (double)Y * HistGram[Y];
    return Sum;
}
//---------------------------------------------------------------------------
double  TCalcThreshold::C(long HistGram[], int Index)
{
    double Sum = 0;
    for (int Y = 0; Y <= Index; Y++)
        Sum += (double)Y * Y * HistGram[Y];
    return Sum;
}
//---------------------------------------------------------------------------
double  TCalcThreshold::D(long HistGram[], int Index)
{
    double Sum = 0;
    for (int Y = 0; Y <= Index; Y++)
        Sum += (double)Y * Y * Y * HistGram[Y];
    return Sum;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetKittlerMinError(long HistGram[])
{
    int X, Y;
    int MinValue, MaxValue;
    int Threshold ;
    int PixelBack, PixelFore;
    double OmegaBack, OmegaFore, MinSigma, Sigma, SigmaBack, SigmaFore;
    for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
    for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
    if (MaxValue == MinValue) return MaxValue; 
    if (MinValue + 1 == MaxValue) return MinValue;
    Threshold = -1;
    MinSigma = 1E+20;
    for (Y = MinValue; Y < MaxValue; Y++)
    {
        PixelBack = 0; PixelFore = 0;
        OmegaBack = 0; OmegaFore = 0;
        for (X = MinValue; X <= Y; X++)
        {
            PixelBack += HistGram[X];
            OmegaBack = OmegaBack + X * HistGram[X];
        }
        for (X = Y + 1; X <= MaxValue; X++)
        {
            PixelFore += HistGram[X];
            OmegaFore = OmegaFore + X * HistGram[X];
        }
        OmegaBack = OmegaBack / PixelBack;
        OmegaFore = OmegaFore / PixelFore;
        SigmaBack = 0; SigmaFore = 0;
        for (X = MinValue; X <= Y; X++) SigmaBack = SigmaBack + (X - OmegaBack) * (X - OmegaBack) * HistGram[X];
        for (X = Y + 1; X <= MaxValue; X++) SigmaFore = SigmaFore + (X - OmegaFore) * (X - OmegaFore) * HistGram[X];
        if (SigmaBack == 0 || SigmaFore == 0)
        {
            if (Threshold == -1)
                Threshold = Y;
        }
        else
        {
            SigmaBack = sqrt(SigmaBack / PixelBack);
            SigmaFore = sqrt(SigmaFore / PixelFore);
            double tmp = SigmaFore / PixelFore ;
            if (tmp<=0) return -1;
            Sigma = 1 + 2 * (PixelBack * log(SigmaBack / PixelBack) + PixelFore * log(tmp));
            if (Sigma < MinSigma)
            {
                MinSigma = Sigma;
                Threshold = Y;
            }
        }
    }
    return Threshold;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetIsoDataThreshold(long HistGram[])
{
    int i, l, toth, totl, h, g = 0;
    for (i = 1; i < 256; i++)
    {
        if (HistGram[i] > 0)
        {
            g = i + 1;
            break;
        }
    }
    while (true)
    {
        l = 0;
        totl = 0;
        for (i = 0; i < g; i++)
        {
            totl = totl + HistGram[i];
            l = l + (HistGram[i] * i);
        }
        h = 0;
        toth = 0;
        for (i = g + 1; i < 256; i++)
        {
            toth += HistGram[i];
            h += (HistGram[i] * i);
        }
        if (totl > 0 && toth > 0)
        {
            l /= totl;
            h /= toth;
            if (g == (int)((l + h) / 2.0))
                break;
        }
        g++;
        if (g > 256 - 2)
        {
            return 0;
        }
    }
    return g;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetShanbhagThreshold(long HistGram[])
{
    int threshold;
    int ih, it;
    int first_bin;
    int last_bin;
    double term;
    double tot_ent;  /* total entropy */
    double min_ent;  /* max entropy */
    double ent_back; /* entropy of the background pixels at a given threshold */
    double ent_obj;  /* entropy of the object pixels at a given threshold */
    double norm_histo[256]; /* normalized histogram */
    double P1[256]; /* cumulative normalized histogram */
    double P2[256];

    int total = 0;
    for (ih = 0; ih < 256; ih++)
        total += HistGram[ih];

    for (ih = 0; ih < 256; ih++)
        norm_histo[ih] = (double)HistGram[ih] / total;

    P1[0] = norm_histo[0];
    P2[0] = 1.0 - P1[0];
    for (ih = 1; ih < 256; ih++)
    {
        P1[ih] = P1[ih - 1] + norm_histo[ih];
        P2[ih] = 1.0 - P1[ih];
    }

    /* Determine the first non-zero bin */
    first_bin = 0;
    for (ih = 0; ih < 256; ih++)
    {
        if (!(abs(P1[ih]) < 2.220446049250313E-16))
        {
            first_bin = ih;
            break;
        }
    }

    /* Determine the last non-zero bin */
    last_bin = 256 - 1;
    for (ih = 256 - 1; ih >= first_bin; ih--)
    {
        if (!(abs(P2[ih]) < 2.220446049250313E-16))
        {
            last_bin = ih;
            break;
        }
    }

    // Calculate the total entropy each gray-level
    // and find the threshold that maximizes it 
    threshold = -1;
    min_ent = (std::numeric_limits<double>::max)();

    for (it = first_bin; it <= last_bin; it++)
    {
        /* Entropy of the background pixels */
        ent_back = 0.0;
        term = 0.5 / P1[it];
        for (ih = 1; ih <= it; ih++)
        { //0+1?
            double tmp = 1.0 - term * P1[ih - 1];
            if (tmp<=0) return -1;
            ent_back -= norm_histo[ih] * log(tmp);
        }
        ent_back *= term;

        /* Entropy of the object pixels */
        ent_obj = 0.0;
        term = 0.5 / P2[it];
        for (ih = it + 1; ih < 256; ih++)
        {
            double tmp = 1.0 - term * P2[ih];
            if (tmp<=0) return -1;
            ent_obj -= norm_histo[ih] * log(tmp);
        }
        ent_obj *= term;

        /* Total entropy */
        tot_ent = abs(ent_back - ent_obj);

        if (tot_ent < min_ent)
        {
            min_ent = tot_ent;
            threshold = it;
        }
    }
    return threshold;
}

//---------------------------------------------------------------------------
long  TCalcThreshold::GetYenThreshold(long HistGram[])
{
    int threshold;
    int ih, it;
    double crit;
    double max_crit;
    double norm_histo[256]; /* normalized histogram */
    double P1[256]; /* cumulative normalized histogram */
    double P1_sq[256];
    double P2_sq[256];

    int total = 0;
    for (ih = 0; ih <256; ih++)
        total += HistGram[ih];

    for (ih = 0; ih < 256; ih++)
        norm_histo[ih] = (double)HistGram[ih] / total;

    P1[0] = norm_histo[0];
    for (ih = 1; ih < 256; ih++)
        P1[ih] = P1[ih - 1] + norm_histo[ih];

    P1_sq[0] = norm_histo[0] * norm_histo[0];
    for (ih = 1; ih < 256; ih++)
        P1_sq[ih] = P1_sq[ih - 1] + norm_histo[ih] * norm_histo[ih];

    P2_sq[256 - 1] = 0.0;
    for (ih = 256 - 2; ih >= 0; ih--)
        P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1];

    /* Find the threshold that maximizes the criterion */
    threshold = -1;
    max_crit =  (std::numeric_limits<double>::min)();
    for (it = 0; it < 256; it++)
	{
        double tmp = P1_sq[it] * P2_sq[it];
        if (tmp<=0) continue;
        double tmp1 = P1[it] * (1.0 - P1[it]);
        if (tmp1<=0) continue;
        crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? log(tmp) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? log(tmp1) : 0.0);
        if (crit > max_crit)
		{
            max_crit = crit;
            threshold = it;
        }
    }
    return threshold;
}
//---------------------------------------------------------------------------
long  TCalcThreshold::GetHuangFuzzyThreshold(long HistGram[])
{
    int X, Y;
    int First, Last;
    int Threshold = -1;
    double BestEntropy = (std::numeric_limits<double>::max)(), Entropy;
    
    for (First = 0; First < 256 && HistGram[First] == 0; First++) ;
    for (Last = 256 - 1; Last > First && HistGram[Last] == 0; Last--) ;
    if (First == Last) return First;
    if (First + 1 == Last) return First;

    std::vector<long> S;
    std::vector<long> W;
    S.resize(Last + 1);
    W.resize(Last + 1);          
    S[0] = HistGram[0];
    for (Y = First > 1 ? First : 1; Y <= Last; Y++)
    {
        S[Y] = S[Y - 1] + HistGram[Y];
        W[Y] = W[Y - 1] + Y * HistGram[Y];
    }

    std::vector<double> Smu;
    Smu.resize(Last + 1 - First);

    for (Y = 1; Y <(int)Smu.size(); Y++)
    {
        double mu = 1 / (1 + (double)Y / (Last - First));

        double tmp = mu;
        if (tmp<=0) return -1;
        double tmp1 = 1-mu;
        if (tmp1<=0) return -1;               

        Smu[Y] = -mu * log(tmp) - (1 - mu) * log(tmp1);
    }
   
    for (Y = First; Y <= Last; Y++)
    {
        Entropy = 0;
        int mu = (int)((double)W[Y] / S[Y]);
        for (X = First; X <= Y; X++)
            Entropy += Smu[abs(X - mu)] * HistGram[X];
        mu = (int)((double)(W[Last] - W[Y]) / (S[Last] - S[Y]));
        for (X = Y + 1; X <= Last; X++)
            Entropy += Smu[abs(X - mu)] * HistGram[X];
        if (BestEntropy > Entropy)
        {
            BestEntropy = Entropy;
            Threshold = Y;
        }
    }
    return Threshold;
}
//---------------------------------------------------------------------------